############################################################
#
#  RNPL program to solve 2 eqns resulting from the sourced 
#  3-D KG wave eqn.
# 
############################################################


#-----------------------------------------------------------
# Definition of memory size (only needed for Fortran)
#-----------------------------------------------------------
system parameter int memsiz := 250000000


#-----------------------------------------------------------
# Definition of parameters and associated default values. 
#-----------------------------------------------------------


#-----------------------------------------------------------
# The following parameters are used in the
# specification of the initial data. 
#-----------------------------------------------------------
parameter float A          := 1.0
parameter float deltax     := 2.0
parameter float deltay     := 2.0
parameter float deltaz     := 2.0
parameter float x0         := 0.0
parameter float y0         := 0.0
parameter float z0         := 0.0
parameter float epsKO      := 0.3 // Kreiss-Oliger dissipation parameter ( 0 <= epsKO < 1)

#-----------------------------------------------------------
# Note that "xmin" and "xmax" are special in that they are 
# also implicitly declared to be parameters via the 
# definition of the coordinate system below.
# Note to self: It turns out you can't declare these in terms
# of the quantities defined above. They must be numbers.
#-----------------------------------------------------------
parameter float xmin       := -100.0
parameter float xmax       :=  100.0
parameter float ymin       := -100.0
parameter float ymax       :=  100.0
parameter float zmin       := -100.0
parameter float zmax       :=  100.0


#-----------------------------------------------------------
# Definition of coordinate system, note that the first 
# coordinate is always assumed to be the time coordinate.
#-----------------------------------------------------------
rect coordinates t, x, y, z


#-----------------------------------------------------------
# Definition of finite-difference grids: [1:Nx] specifies
# the index range, {xmin:xmax} the coordinate range.
#-----------------------------------------------------------
uniform rect[x,y,z] grid g1 [1:Nx][1:Ny][1:Nz] {xmin:xmax}{ymin:ymax}{zmin:zmax}

#-----------------------------------------------------------
# Definition of grid functions: since a Crank-Nicholson 
# scheme is being used, the grid functions are defined at 
# "temporal offsets" 0 and 1, corresponding to "current" and 
# "advanced" time levels.  The directive {out_gf = 1} enables
# default output of the grid function (output can be disabled
# via {out_gf = 0}, by omitting the directive completely, or 
# via modification of the file .rnpl.attributes prior to 
# program invocation).
#-----------------------------------------------------------
float phi  on g1 at 0,1 {out_gf = 1}
float mu   on g1 at 0,1 {out_gf = 1}

#Right hand sides of the residual equations
float RHSphi on g1 at 0,1 {out_gf = 0}
float RHSmu on g1 at 0,1 {out_gf = 0}


#-----------------------------------------------------------
# FINITE DIFFERENCE OPERATOR DEFINITIONS
#-----------------------------------------------------------


#-----------------------------------------------------------
# Forward in time operator (for sanity check use)
#-----------------------------------------------------------
operator FORW(f,t) := <1>f[0][0][0]

#-----------------------------------------------------------
# Crank Nicholson time derivative operator (first forward 
# difference)
#-----------------------------------------------------------
operator DCN(f,t) := (<1>f[0][0][0] - <0>f[0][0][0]) / dt    



#-----------------------------------------------------------
# Kreiss-Oliger dissipation operator
#-----------------------------------------------------------
operator DISS_KO(f,t) := - (epsKO/(16*dt))* (<0>f[-2][0][0] - 4*<0>f[-1][0][0] + 6*<0>f[0][0][0] - 4*<0>f[1][0][0] + <0>f[2][0][0] + 
                                             <0>f[0][-2][0] - 4*<0>f[0][-1][0] + 6*<0>f[0][0][0] - 4*<0>f[0][1][0] + <0>f[0][2][0] +
                                             <0>f[0][0][-2] - 4*<0>f[0][0][-1] + 6*<0>f[0][0][0] - 4*<0>f[0][0][1] + <0>f[0][0][2]
                                           )



#-----------------------------------------------------------
# DIFFERENCE EQUATION DEFINITIONS
# This section has been disabled, as the equations have been
# written straight in the updates.f and residual.f files.
#-----------------------------------------------------------

evaluate residual phi 
{ 
     [1:1][1:Ny][1:Nz]        := FORW(phi,t)  = 0;
     [1:Nx][1:1][1:Nz]        := FORW(phi,t)  = 0;
     [1:Nx][1:Ny][1:1]        := FORW(phi,t)  = 0;

     [2:2][2:Ny-1][2:Nz-1]    := DCN(phi,t)   = RHSphi;
     [2:Nx-1][2:2][2:Nz-1]    := DCN(phi,t)   = RHSphi;    
     [2:Nx-1][2:Ny-1][2:2]    := DCN(phi,t)   = RHSphi;  

     [3:Nx-2][3:Ny-2][3:Nz-2] := DCN(phi,t)   = DISS_KO(phi,t) +RHSphi;

     [Nx-1:Nx-1][2:Ny-1][2:Nz-1]    := DCN(phi,t)   = RHSphi;
     [2:Nx-1][Ny-1:Ny-1][2:Nz-1]    := DCN(phi,t)   = RHSphi;     
     [2:Nx-1][2:Ny-1][Nz-1:Nz-1]    := DCN(phi,t)   = RHSphi;  

     [Nx:Nx][1:Ny][1:Nz]      := FORW(phi,t)  = 0;
     [1:Nx][Ny:Ny][1:Nz]      := FORW(phi,t)  = 0;
     [1:Nx][1:Ny][Nz:Nz]      := FORW(phi,t)  = 0;
}

evaluate residual mu 
{ 
     [1:1][1:Ny][1:Nz]        := FORW(mu,t)  = 0;
     [1:Nx][1:1][1:Nz]        := FORW(mu,t)  = 0;
     [1:Nx][1:Ny][1:1]        := FORW(mu,t)  = 0;

     [2:2][2:Ny-1][2:Nz-1]    := DCN(mu,t)  = RHSmu;
     [2:Nx-1][2:2][2:Nz-1]    :=  DCN(mu,t) = RHSmu;   
     [2:Nx-1][2:Ny-1][2:2]    := DCN(mu,t)  = RHSmu;

     [3:Nx-2][3:Ny-2][3:Nz-2] :=  DCN(mu,t) = DISS_KO(mu,t) +RHSmu;

     [Nx-1:Nx-1][2:Ny-1][2:Nz-1]    := DCN(mu,t)  = RHSmu;
     [2:Nx-1][Ny-1:Ny-1][2:Nz-1]    := DCN(mu,t)  = RHSmu;    
     [2:Nx-1][2:Ny-1][Nz-1:Nz-1]    := DCN(mu,t)  = RHSmu;

     [Nx:Nx][1:Ny][1:Nz]      := FORW(mu,t)  = 0;
     [1:Nx][Ny:Ny][1:Nz]      := FORW(mu,t)  = 0;
     [1:Nx][1:Ny][Nz:Nz]      := FORW(mu,t)  = 0;
}


#-----------------------------------------------------------
# INITIALIZATION STATEMENTS 
#-----------------------------------------------------------

initialize phi 
{
     [1:Nx][1:Ny][1:Nz] := A* exp( -(x-x0)^2 / deltax^2 -(y-y0)^2 / deltay^2 -(z-z0)^2 / deltaz^2 )
}

initialize mu 
{
    [1:Nx][1:Ny][1:Nz] := - A* exp( -(x-x0)^2 / deltax^2 -(y-y0)^2 / deltay^2 -(z-z0)^2 / deltax^2 )
}

initialize RHSphi
{
   [1:Nx][1:Ny][1:Nz] := 0;
}

initialize RHSmu
{
   [1:Nx][1:Ny][1:Nz] := 0;
}


auto initialize phi, mu, RHSphi, RHSmu


#-----------------------------------------------------------
# Definition of type of time stepping algorithm.  The use
# of "iterative" here, combined with the "auto update ..." 
# statement below, results in a scheme whereby the 
# residuals defined above are iteratively relaxed using 
# a point-wise Newton-Gauss-Seidel technique, until the 
# residual norms are below a certain threshold.
#-----------------------------------------------------------
looper iterative

#-----------------------------------------------------------
# The following statement directs RNPL to automatically 
# generate code to update the grid functions using the 
# residual definitions.
#-----------------------------------------------------------

auto update phi, mu

rhs.inc update_rhs updates RHSphi, RHSmu
   header RHSphi, RHSmu, phi, mu, t, x, y, z, dx, dy, dz, dt
